home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / eigenvec.pro < prev    next >
Text File  |  1997-07-08  |  9KB  |  207 lines
  1. ;$Id: eigenvec.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       EIGENVEC
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the eigenvectors of an N by N real, non-
  11. ;       symmetric array using inverse subspace iteration. The result is 
  12. ;       a complex array with a column dimension equal to N and a row 
  13. ;       dimension equal to the number of eigenvalues.
  14. ;
  15. ; CATEGORY:
  16. ;       Linear Algebra / Eigensystems
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = Eigenvec(A, Eval)
  20. ;
  21. ; INPUTS:
  22. ;       A:    An N by N nonsymmetric array of type float or double.
  23. ;
  24. ;    EVAL:    An N-element complex vector of eigenvalues.
  25. ;
  26. ; KEYWORD PARAMETERS:
  27. ;       DOUBLE:  If set to a non-zero value, computations are done in
  28. ;                double precision arithmetic.
  29. ;
  30. ;        ITMAX:  The number of iterations performed in the computation
  31. ;                of each eigenvector. The default value is 4.
  32. ;
  33. ;     RESIDUAL:  Use this keyword to specify a named variable which returns
  34. ;                the residuals for each eigenvalue/eigenvector(lambda/x) pair.
  35. ;                The residual is based on the definition Ax - (lambda)x = 0
  36. ;                and is an array of the same size and type as RESULT. The rows
  37. ;                this array correspond to the residuals for each eigenvalue/
  38. ;                eigenvector pair. This keyword must be initialized to a non-
  39. ;                zero value before calling EIGENVEC if the residuals are 
  40. ;                desired.
  41. ;
  42. ; EXAMPLE:
  43. ;       Define an N by N real, nonsymmetric array.
  44. ;         a = [[1.0, -2.0, -4.0,  1.0], $
  45. ;              [0.0, -2.0,  3.0,  4.0], $
  46. ;              [2.0, -6.0, -1.0,  4.0], $
  47. ;              [3.0, -3.0,  1.0, -2.0]]
  48. ;
  49. ;       Compute the eigenvalues of A using double-precision complex arithmetic.
  50. ;         eval = HQR(ELMHES(a), /double)
  51. ;
  52. ;       Print the eigenvalues. The correct solution should be:
  53. ;       (0.26366259, -6.1925899), (0.26366259, 6.1925899), $
  54. ;       (-4.9384492,  0.0000000), (0.41112397, 0.0000000)
  55. ;         print, eval
  56. ;
  57. ;       Compute the eigenvectors of A. The eigenvectors are returned in the 
  58. ;       rows of EVEC. The residual keyword must be initialized as a nonzero
  59. ;       value prior to calling EIGENVEC.
  60. ;         residual = 1
  61. ;         result = EIGENVEC(a, eval, residual = residual)
  62. ;
  63. ;       Print the eigenvectors.
  64. ;         print, evec(*,0), evec(*,1), evec(*,2), evec(*,3)
  65. ;
  66. ;       The accuracy of each eigenvalue/eigenvector (lamda/x) 
  67. ;       pair may be checked by printing the residual array. This array is the
  68. ;       same size and type as RESULT and returns the residuals as its rows.
  69. ;       The residual is based on the mathematical definition of an eigenvector,
  70. ;       Ax - (lambda)x = 0. All residual values should be floating-point zeros.
  71. ;         print, residual
  72. ;
  73. ; PROCEDURE:
  74. ;       EIGENVEC computes the set of eigenvectors that correspond to a given 
  75. ;       set of eigenvalues using Inverse Subspace Iteration. The eigenvectors 
  76. ;       are computed up to a scale factor and are of Euclidean length. The
  77. ;       existence and uniqueness of eigenvectors are not guaranteed.
  78. ;
  79. ; MODIFICATION HISTORY:
  80. ;       Written by:  GGS, RSI, December 1994
  81. ;       Modified:    GGS, RSI, April 1996
  82. ;                    Modified keyword checking and use of double precision. 
  83. ;-
  84.  
  85. FUNCTION EigenVec, A, Eval, Double = Double, ItMax = ItMax, $
  86.                                              Residual = Residual
  87.  
  88.   ON_ERROR, 2  ;Return to caller if error occurs.
  89.  
  90.   if N_PARAMS() ne 2 then $
  91.     MESSAGE, "Incorrect number of input arguments."
  92.     
  93.   TypeA = SIZE(A)
  94.   TypeEval = SIZE(Eval)
  95.  
  96.   if TypeA[1] ne TypeA[2] then $
  97.     MESSAGE, "Input array must be square."
  98.  
  99.   if TypeA[3] ne 4 and TypeA[3] ne 5 then $
  100.     MESSAGE, "Input array must be float or double."
  101.  
  102.   if TypeEval[TypeEval[0]+1] ne 6 and TypeEval[TypeEval[0]+1] ne 9 then $
  103.     MESSAGE, "Eigenvalues must be complex or double-complex."
  104.  
  105.   Enum = TypeEval[TypeEval[0]+2] ;Number of eigenvalues.
  106.   if TypeA[2] ne Enum then $
  107.     MESSAGE, "Input array and eigenvalues are of incompatible size."
  108.  
  109.   ;If the DOUBLE keyword is not set then the internal precision and
  110.   ;result are determined by the type of input.
  111.   if N_ELEMENTS(Double) eq 0 then $
  112.     Double = (TypeA[TypeA[0]+1] eq 5 or TypeEval[TypeEval[0]+1] eq 9)
  113.  
  114.   if N_ELEMENTS(ItMax) eq 0 then ItMax = 4
  115.  
  116.   Diag = LINDGEN(TypeA[1]) * (TypeA[1]+1) ;Diagonal indices.
  117.  
  118.   ;Double Precision.
  119.   if Double ne 0 then begin
  120.     Evec = DCOMPLEXARR(TypeA[1], Enum) ;Eigenvector storage array with number
  121.                                     ;of rows equal to number of eigenvalues.
  122.     for k = 0, Enum - 1 do begin
  123.       Alud = A  ;Create a copy of the array for next eigenvalue computation.
  124.       if IMAGINARY(Eval[k]) ne 0 then begin ;Complex eigenvalue.
  125.         Alud = DCOMPLEX(Alud)
  126.         Alud[Diag] = Alud[Diag] - Eval[k]
  127.         ;Avoid intermediate variables. re = DOUBLE(Alud) im = IMAGINARY(Alud)
  128.         Comp = [[DOUBLE(Alud), -IMAGINARY(Alud)], $
  129.                 [IMAGINARY(Alud), DOUBLE(Alud)]]
  130.         ;Initial eigenvector.
  131.         B = REPLICATE(1.0d, 2*TypeA[1]) / SQRT(2.0d * TypeA[1])
  132.         LUDC, Comp, Index, DOUBLE = DOUBLE
  133.         it = 0
  134.         while it lt ItMax do begin ;Iteratively compute the eigenvector.
  135.           X = LUSOL(Comp, Index, B, DOUBLE = DOUBLE)
  136.           B = X / SQRT(TOTAL(X^2, 1, DOUBLE = DOUBLE)) ;Normalize eigenvector.
  137.           it = it + 1
  138.         endwhile
  139.         ;Row vector storage.
  140.         Evec[*, k] = DCOMPLEX(B[0:TypeA[1]-1], B[TypeA[1]:*])
  141.       endif else begin ;Real eigenvalue
  142.         Alud[Diag] = Alud[Diag] - DOUBLE(Eval[k])
  143.         B = REPLICATE(1.0d, TypeA[1]) / SQRT(TypeA[1]+0.0d)
  144.         LUDC, Alud, Index, DOUBLE = DOUBLE
  145.         it = 0
  146.         while it lt ItMax do begin
  147.           X = LUSOL(Alud, Index, B, DOUBLE = DOUBLE)
  148.           B = X / SQRT(TOTAL(X^2, 1, DOUBLE = DOUBLE)) ;Normalize eigenvector.
  149.           it = it + 1
  150.         endwhile
  151.         Evec[*, k] = DCOMPLEX(B, 0.0d0) ;Row vector storage.
  152.       endelse
  153.     endfor
  154.     if keyword_set(Residual) then begin ;Compute eigenvalue/vector residuals.
  155.       ;Because RESIDUAL is a keyword that envokes functionality and returns
  156.       ;a named variable, it must be initialized before calling EIGENVEC().
  157.       Residual = DCOMPLEXARR(TypeA[1], Enum) ;Dimensioned the same as Evec.
  158.         for k = 0, Enum - 1 do $
  159.           Residual[*,k] = (A##Evec[*,k]) - (Eval[k] * Evec[*,k])
  160.     endif
  161.   endif else begin ;Single Precision.
  162.     Evec = COMPLEXARR(TypeA[1], Enum) ;Eigenvector storage array.
  163.     for k = 0, Enum - 1 do begin
  164.       Alud = A  ;Create a copy of the array for next eigenvalue computation.
  165.       if IMAGINARY(Eval[k]) ne 0 then begin ;Complex eigenvalue.
  166.         Alud = COMPLEX(Alud)
  167.         Alud[Diag] = Alud[Diag] - Eval[k]
  168.         ;Avoid intermediate variables. re = FLOAT(Alud) im = IMAGINARY(Alud)
  169.         Comp = [[FLOAT(Alud), -IMAGINARY(Alud)], $
  170.                 [IMAGINARY(Alud), FLOAT(Alud)]]
  171.         ;Initial eigenvector.
  172.         B = REPLICATE(1.0, 2*TypeA[1]) / SQRT(2.0 * TypeA[1])
  173.         LUDC, Comp, Index, DOUBLE = DOUBLE
  174.         it = 0
  175.         while it lt ItMax do begin ;Iteratively compute the eigenvector. 
  176.           X = LUSOL(Comp, Index, B, DOUBLE = DOUBLE)
  177.           B = X / SQRT(TOTAL(X^2, 1)) ;Normalize eigenvector.
  178.           it = it + 1
  179.         endwhile
  180.         ;Row vector storage.
  181.         Evec[*, k] = COMPLEX(B[0:TypeA[1]-1], B[TypeA[1]:*])
  182.       endif else begin ;Real eigenvalue 
  183.         Alud[Diag] = Alud[Diag] - FLOAT(Eval[k])
  184.         B = REPLICATE(1.0, TypeA[1]) / SQRT(TypeA[1])
  185.         LUDC, Alud, Index, DOUBLE = DOUBLE
  186.         it = 0
  187.         while it lt ItMax do begin
  188.           X = LUSOL(Alud, Index, B, DOUBLE = DOUBLE)
  189.           B = X / SQRT(TOTAL(X^2, 1))  ;Normalize eigenvector.
  190.           it = it + 1
  191.         endwhile
  192.         Evec[*, k] = COMPLEX(B, 0.0) ;Row vector storage.
  193.       endelse
  194.     endfor
  195.     if keyword_set(Residual) then begin ;Compute eigenvalue/vector residuals.
  196.       ;Because RESIDUAL is a keyword that envokes functionality and returns
  197.       ;a named variable, it must be initialized before calling EIGENVEC().
  198.       Residual = COMPLEXARR(TypeA[1], Enum) ;Dimensioned the same as Evec.
  199.         for k = 0, Enum - 1 do $
  200.           Residual[*,k] = (A##Evec[*,k]) - (Eval[k] * Evec[*,k])
  201.     endif
  202.   endelse
  203.   
  204.   if Double eq 0 then RETURN, COMPLEX(Evec) else RETURN, Evec
  205.  
  206. END
  207.